{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ab8f1e6c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "IPython.notebook.set_autosave_interval(120000)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Autosaving every 120 seconds\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "%autosave 120"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d3a38e45",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from functools import reduce\n",
    "from Matrix_Solutions import Matrix_Solutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a4b93330",
   "metadata": {},
   "outputs": [],
   "source": [
    "MS=Matrix_Solutions()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d707ff35",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2., -1.,  1.],\n",
       "       [-1., -2.,  3.],\n",
       "       [ 1.,  3.,  1.]], dtype=float32)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A=np.array([[2,-1,1],[-1,-2,3],[1,3,1]],dtype=np.float32).reshape((3,3))\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "87f720ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "def domain_row_transform(arr,eps=1e-6,Test=False):\n",
    "    assert len(arr.shape)==2\n",
    "    A = np.copy(arr)\n",
    "    M, N = A.shape\n",
    "    for i in range(M):\n",
    "        for j in range(i, N):\n",
    "            if Test:\n",
    "                print('During process in row:{}, col:{}'.format(i, j))\n",
    "            if sum(abs(A[i:, j])) > eps:\n",
    "                zero_index = []\n",
    "                non_zero_index = []\n",
    "                for k in range(i, M):\n",
    "                    if abs(A[k, j]) > eps:\n",
    "                        non_zero_index.append(k)\n",
    "                    else:\n",
    "                        zero_index.append(k)\n",
    "\n",
    "                non_zero_index = sorted(\n",
    "                    non_zero_index, key=lambda x: abs(A[x, j]))\n",
    "                sort_idnex = non_zero_index+zero_index\n",
    "\n",
    "                if Test:\n",
    "                    print('Sorting index for {} cols:\\n'.format(i), sort_idnex)\n",
    "\n",
    "                if sort_idnex[0] != i:\n",
    "                    A[[sort_idnex[0], i], :] = A[[i, sort_idnex[0]], :]\n",
    "                if Test:\n",
    "                    print('After resorting\\n', A)\n",
    "    return A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9e815f6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ddbbf7dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Doolittle_Decomposition(target, domain_selection=False, test=False):\n",
    "    \"\"\"Doolittle_Decomposition equal to LU_Decomposition\n",
    "\n",
    "    Args:\n",
    "        target ([np.darray]): [input matrix]\n",
    "        domain_selection (bool, optional): [operate domain selection on target]. Defaults to True.\n",
    "        test (bool, optional): [show testing information]. Defaults to False.\n",
    "\n",
    "    Returns:\n",
    "        [resort input, L, U] if domain_selection=True else [L,U]\n",
    "    \"\"\"\n",
    "    assert len(target.shape) == 2\n",
    "    arr = np.copy(target)\n",
    "    M, N = arr.shape\n",
    "    L = np.zeros_like(arr, dtype=arr.dtype)\n",
    "    U = np.zeros_like(arr, dtype=arr.dtype)\n",
    "    if domain_selection:\n",
    "        arr = domain_row_transform(arr)\n",
    "    for i in range(M):\n",
    "        if test:\n",
    "            print(i)\n",
    "            print(L, U, sep='\\n')\n",
    "\n",
    "        for j in range(i, N):\n",
    "            if test:\n",
    "                print('U', i, j, arr[i, j], L[i, :i], U[:i, j])\n",
    "            U[i, j] = arr[i, j]-L[i, :i]@U[:i, j]\n",
    "            if i == j:\n",
    "                L[i, j] = 1\n",
    "            else:\n",
    "                if test:\n",
    "                    print('L', j, i, arr[j, i],\n",
    "                          L[j, :i], U[:i, i], U[i, i])\n",
    "                L[j, i] = (arr[j, i]-(L[j, :i]@U[:i, i]))/U[i, i]\n",
    "    if domain_selection:\n",
    "        print('return resort input matrix, L, U')\n",
    "        return arr, L, U\n",
    "    else:\n",
    "        return L, U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5f51865d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.50948856, -1.61033161, -0.20671825,  0.4286645 , -0.2572608 ],\n",
       "       [-1.28615783,  0.87991375,  1.71080037, -1.61832726,  1.52242031],\n",
       "       [-0.98092031, -0.32273876,  2.20229521,  0.853808  , -2.15168542],\n",
       "       [-2.02003405,  0.05052931, -1.01116616,  0.26437957, -1.65870163],\n",
       "       [-0.83327297,  0.32495716, -1.02878403, -0.96607898, -0.59284746]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A=np.random.randn(5,5)\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 266,
   "id": "74300a3e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.165  0.618 -0.885 -0.515 -1.011]\n",
      " [-0.774  2.898 -0.716  0.118  1.437]\n",
      " [-0.905  0.983 -1.104 -0.89   0.997]\n",
      " [-1.395  0.064  0.171 -0.275  0.52 ]\n",
      " [ 1.446  0.897 -1.    -2.323  0.294]]\n",
      "[[    1.        0.        0.        0.        0.   ]\n",
      " [    4.688     1.        0.        0.        0.   ]\n",
      " [    5.484 -3863.106     1.        0.        0.   ]\n",
      " [    8.449 -8279.836     2.143     1.        0.   ]\n",
      " [   -8.763 10132.862    -2.623   -11.842     1.   ]]\n",
      "[[   -0.165     0.618    -0.885    -0.515    -1.011]\n",
      " [    0.        0.001     3.434     2.532     6.176]\n",
      " [    0.        0.    13269.089  9782.795 23865.142]\n",
      " [    0.        0.        0.        0.216    -4.26 ]\n",
      " [    0.        0.        0.        0.      -43.798]]\n"
     ]
    }
   ],
   "source": [
    "TA,L1,U1=MS.LDV_2D(A,mode=\"LU\")\n",
    "print(TA,L1,U1,sep='\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "75aa45e5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.          0.          0.          0.          0.        ]\n",
      " [ 0.85204874  1.          0.          0.          0.        ]\n",
      " [ 0.6498362   0.32136532  1.          0.          0.        ]\n",
      " [ 1.33822415  0.97936016 -1.49258468  1.          0.        ]\n",
      " [ 0.55202337  0.53903227 -1.11649086  0.35443372  1.        ]]\n",
      "[[-1.50948856 -1.61033161 -0.20671825  0.4286645  -0.2572608 ]\n",
      " [ 0.          2.25199476  1.8869344  -1.9835703   1.74161905]\n",
      " [ 0.          0.          1.73023294  1.21269699 -2.544204  ]\n",
      " [ 0.          0.          0.          3.44341308 -6.81754125]\n",
      " [ 0.          0.          0.          0.         -1.81383636]]\n"
     ]
    }
   ],
   "source": [
    "L,U=Doolittle_Decomposition(A,test=False)\n",
    "print(L,U,sep='\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "8dfa665d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.50948856, -1.61033161, -0.20671825,  0.4286645 , -0.2572608 ],\n",
       "       [-1.28615783,  0.87991375,  1.71080037, -1.61832726,  1.52242031],\n",
       "       [-0.98092031, -0.32273876,  2.20229521,  0.853808  , -2.15168542],\n",
       "       [-2.02003405,  0.05052931, -1.01116616,  0.26437957, -1.65870163],\n",
       "       [-0.83327297,  0.32495716, -1.02878403, -0.96607898, -0.59284746]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L@U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 274,
   "id": "f8381cdc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0., -0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.]])"
      ]
     },
     "execution_count": 274,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L@U-L1@U1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 253,
   "id": "dd751085",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.099,  1.973,  1.182, -0.766,  0.645],\n",
       "       [-0.217, -0.395,  1.096,  0.408,  1.363],\n",
       "       [ 0.267, -1.086,  0.515, -0.484, -0.904],\n",
       "       [ 1.049,  0.281, -0.146, -0.868,  0.975],\n",
       "       [-3.133,  0.649,  1.842, -1.076,  0.508]])"
      ]
     },
     "execution_count": 253,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 237,
   "id": "4f4bea3e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.391,  3.218,  0.912, -2.658,  1.12 ],\n",
       "       [ 0.942, -1.41 ,  0.645, -0.643,  0.385],\n",
       "       [ 1.936, -0.036, -1.194,  1.028,  3.794],\n",
       "       [-0.863, -1.303, -0.303,  2.748, -2.38 ],\n",
       "       [-1.624, -0.469, -0.061, -0.475, -2.919]])"
      ]
     },
     "execution_count": 237,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L@U-nA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "bf1d16b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_SDD_M(arr,eps=1e-4):\n",
    "    \"\"\"check_SDDM 判断输入矩阵是否满足严格对角占优条件\n",
    "    Args:\n",
    "        arr ([np.darray]): [input matrix]\n",
    "        eps ([type], optional): [threshold]. Defaults to 1e-4.\n",
    "\n",
    "    Returns:\n",
    "        [bool]: [whether or not]\n",
    "    \"\"\"\n",
    "    assert arr.shape[0]==arr.shape[1]\n",
    "    return check_SDD_vec(np.diagonal(arr,offset=-1),np.diagonal(arr,offset=0),np.diagonal(arr,offset=1),eps=eps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "a24eb255",
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_SDD_vec(a,b,c,eps=1e-4,test=False):\n",
    "    assert len(a)==len(c) and len(a)+1 == len(b)\n",
    "    N=len(b)\n",
    "    a_=np.concatenate(([eps],a),axis=0)\n",
    "    c_=np.concatenate((c,[eps]),axis=0)\n",
    "    if test:\n",
    "        print(a_,c_)\n",
    "    sum_a_c=np.abs(a_)+np.abs(c_)\n",
    "    b_=np.abs(b)\n",
    "    for i in range(N):\n",
    "        if b_[i]<sum_a_c[i]:\n",
    "            return False\n",
    "    return True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "865c4877",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[3., 1., 0., 0.],\n",
       "       [1., 2., 1., 0.],\n",
       "       [0., 1., 3., 2.],\n",
       "       [0., 0., 3., 4.]], dtype=float32)"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A=np.array([[3,1,0,0],[1,2,1,0],[0,1,3,2],[0,0,3,4]],dtype=np.float32).reshape((4,4))\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "3eac3eca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "check_SDD_vec(np.diagonal(A,offset=-1),np.diagonal(A,offset=0),np.diagonal(A,offset=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "c1a7fa2a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "check_SDD_M(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "296b240f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.35812746]\n",
      " [ 0.44621921]\n",
      " [-0.45669147]\n",
      " [ 1.85002579]]\n",
      "[[-0.62816316]\n",
      " [ 0.07761949]\n",
      " [ 2.77619638]\n",
      " [ 6.03002876]]\n"
     ]
    }
   ],
   "source": [
    "x=np.random.randn(4,1)\n",
    "print(x)\n",
    "temp=A@x\n",
    "print(temp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "47805336",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "check_SDDM(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "a1e6f067",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Thomas_Linear_Equation_Solve(a,b,c,d,test=False):\n",
    "    \"\"\"Thomas_Linear_Equation_Solve \n",
    "\n",
    "    Args:\n",
    "        a ([list]): [lower-secondary diagonal]\n",
    "        b ([list]): [diagonal]\n",
    "        c ([list]): [upper-secondary diagonal]\n",
    "        d ([list]): [b of Ax=b]\n",
    "        test (bool, optional): [show testing information]. Defaults to False.\n",
    "\n",
    "    Returns:\n",
    "        x [list]: [solution of Ax=b]\n",
    "        \n",
    "    Author: Junno\n",
    "    \n",
    "    Date: 2022-04-28\n",
    "    \"\"\"\n",
    "    assert len(a)==len(c) and len(a)+1 == len(b)\n",
    "    if not check_SDD_vec(a,b,c):\n",
    "        print('The answer may be wrong because the input are not meet SDD condition')\n",
    "    N=len(b)\n",
    "    if test:\n",
    "        print(a,b,c,d)\n",
    "    beta=np.zeros_like(b)\n",
    "    x=np.zeros_like(b)\n",
    "    y=np.zeros_like(b)\n",
    "    beta[0]=c[0]/b[0]\n",
    "    y[0]=d[0]/b[0]\n",
    "    for i in range(1,N):\n",
    "        if i<N-1:\n",
    "            beta[i]=c[i]/(b[i]-a[i-1]*beta[i-1])\n",
    "        y[i]=(d[i]-a[i-1]*y[i-1])/(b[i]-a[i-1]*beta[i-1])\n",
    "    for i in range(N-1,-1,-1):\n",
    "        if i==N-1:\n",
    "            x[i]=y[i]\n",
    "        else:\n",
    "            x[i]=y[i]-beta[i]*x[i+1]\n",
    "    \n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "d00d3168",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 1 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-0.35812744,  0.44621915, -0.45669138,  1.8500258 ], dtype=float32)"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "Thomas_Linear_Equation_Solve(np.diagonal(A,offset=-1),np.diagonal(A,offset=0),np.diagonal(A,offset=1),temp,test=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "69eb30b4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.35812746],\n",
       "       [ 0.44621921],\n",
       "       [-0.45669147],\n",
       "       [ 1.85002579]])"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "id": "536346ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Matrix_Solutions import Matrix_Solutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "id": "c049edb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "MS=Matrix_Solutions()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "id": "b8d53765",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 1.04 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 0.424],\n",
       "       [-0.909],\n",
       "       [-0.524],\n",
       "       [ 0.451]])"
      ]
     },
     "execution_count": 163,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "MS.Solve_Ax_b_Equation(A,temp,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "id": "721b02f1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.919],\n",
       "       [-1.581],\n",
       "       [ 0.23 ],\n",
       "       [-0.546]])"
      ]
     },
     "execution_count": 164,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "temp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9137993",
   "metadata": {},
   "source": [
    "# Jacobi迭代法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 221,
   "id": "bccefbe3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.862,  0.196, -1.52 , -0.654],\n",
       "       [-0.107, -0.597,  0.842, -0.861],\n",
       "       [ 0.007,  1.341, -1.022, -2.181],\n",
       "       [ 0.981,  0.497, -0.423,  0.547]])"
      ]
     },
     "execution_count": 221,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# A=np.array([[3,1,0,0],[1,2,1,0],[0,1,3,2],[0,0,3,4]],dtype=np.float32).reshape((4,4))\n",
    "A=np.random.randn(4,4)\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "52a6470f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.178]\n",
      " [-0.471]\n",
      " [-0.099]\n",
      " [-0.02 ]]\n",
      "[[-0.04 ]\n",
      " [-1.162]\n",
      " [-0.622]\n",
      " [-0.332]]\n"
     ]
    }
   ],
   "source": [
    "x=np.random.randn(A.shape[0],1)\n",
    "print(x)\n",
    "b=A@x\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "ec6c02c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "U=np.triu(A,k=1)\n",
    "L=np.tril(A,k=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "id": "56811b8e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.905 0.539 0.099 0.023]\n"
     ]
    }
   ],
   "source": [
    "D=np.diag(np.diag(A))\n",
    "B=-np.linalg.inv(D)@(L+U)\n",
    "u,s,vh=np.linalg.svd(B)\n",
    "print(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "333373c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "max_e=1\n",
    "while max_e>=1:\n",
    "    A=np.random.randn(4,4)\n",
    "    U=np.triu(A,k=1)\n",
    "    L=np.tril(A,k=-1)\n",
    "    D=np.diag(np.diag(A))\n",
    "    B=-np.linalg.inv(D)@(L+U)\n",
    "    u,s,vh=np.linalg.svd(B)\n",
    "    max_e=np.max(np.diag(s))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "id": "27da71fe",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.758,  0.28 ,  0.383,  0.259],\n",
       "       [-0.165,  2.509, -0.542,  0.271],\n",
       "       [ 0.762,  1.151,  1.999,  0.842],\n",
       "       [-0.724,  0.515, -0.081, -1.539]])"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "6fed2760",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.35 ,  0.209, -0.218, -0.123],\n",
       "       [ 0.179,  1.616, -0.104, -0.497],\n",
       "       [ 0.333,  0.278,  2.417, -0.226],\n",
       "       [ 1.128,  0.346,  0.104,  1.237]], dtype=float32)"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A=np.array([[-1.35 ,  0.209, -0.218, -0.123],\n",
    "       [ 0.179,  1.616, -0.104, -0.497],\n",
    "       [ 0.333,  0.278,  2.417, -0.226],\n",
    "       [ 1.128,  0.346,  0.104,  1.237]],dtype=np.float32).reshape((4,4))\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "6c8e4e89",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'MS' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m~\\AppData\\Local\\Temp/ipykernel_6332/924765697.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mu\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0ms\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mMS\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msvd\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mB\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mtest\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ms\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'MS' is not defined"
     ]
    }
   ],
   "source": [
    "u,s,v=MS.svd(B,test=False)\n",
    "np.diag(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "id": "854fec45",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Serial_iteration(A, b, N=0, eps=1e-4, x_0=None, iter_way='Gauss_Seidel', norm='infty', own_method=False, test=False):\n",
    "    \"\"\"Serial_iteration:Jacobi and Gauss_Seidel\n",
    "\n",
    "    Args:\n",
    "        A ([np.darray]): [A]\n",
    "        b ([np.darray]): [b]\n",
    "        N (int, optional): [iteration time, it will give a suitable one when N is 0]. Defaults to 0.\n",
    "        eps ([float], optional): [error threshold]. Defaults to 1e-4.\n",
    "        x_0 ([np.darray], optional): [initial x for iteration]. Defaults to None.\n",
    "        norm (str, optional): [iteration method,'Gauss_Seidel','Jacobi']. Defaults to 'Gauss_Seidel'.\n",
    "        norm (str, optional): [norm method to calculate fitting error,'infty','L1']. Defaults to 'infty'.\n",
    "        own_method (bool, optional): [use own method to calculate inv or norm, etc]. Defaults to False.\n",
    "        test (bool, optional): [show testing information]. Defaults to False.\n",
    "\n",
    "    Raises:\n",
    "        ValueError: [The spetrum of Iteration Matrix must be smaller than 1]\n",
    "\n",
    "    Returns:\n",
    "        [x]: [solution of Ax=b]\n",
    "    \n",
    "    Author: Junno\n",
    "        \n",
    "    Date: 2022-04-29\n",
    "    \"\"\" \n",
    "    assert A.shape[0] == b.shape[0]\n",
    "    # use np method or your own\n",
    "    Inv = np.linalg.inv if not own_method else Matrix_Solutions.Primary_Row_Transformation_Inv\n",
    "    SVD = np.linalg.svd if not own_method else Matrix_Solutions.svd\n",
    "    # substract D,L,U from A\n",
    "    D, L, U = np.diag(np.diag(A)), np.tril(A, k=-1), np.triu(A, k=1)\n",
    "    # calculate B and f\n",
    "    if iter_way == 'Jacobi':\n",
    "        B = -Inv(D)@(L+U)\n",
    "        f = Inv(D)@b\n",
    "    elif iter_way == 'Gauss_Seidel':\n",
    "        B = -Inv(D+L)@U\n",
    "        f = Inv(D+L)@b\n",
    "    # initial x for iteration\n",
    "    if x_0 == None:\n",
    "        x_0 = np.zeros_like(b)\n",
    "    \n",
    "    # get eig of B\n",
    "    u, s, v = SVD(B)\n",
    "    # spectrum of B\n",
    "    eig_B = np.diag(s)\n",
    "    max_eig_B = np.max(eig_B)\n",
    "    # check convergence condition\n",
    "    if max_eig_B > 1:\n",
    "        print('The spetrum of Iteration Matrix must be smaller than 1')\n",
    "        raise ValueError\n",
    "        \n",
    "    # matrix norm\n",
    "    if norm == 'infty':\n",
    "        F = Matrix_Solutions.m_infinite_norm\n",
    "    elif norm == 'L1':\n",
    "        F = Matrix_Solutions.m_L1_norm\n",
    "\n",
    "    q = F(B)\n",
    "    if test:\n",
    "        print('q =', q)\n",
    "    \n",
    "    # iteration times\n",
    "    N = int(np.ceil(np.log10(1/eps)*np.log(10) /\n",
    "            (-np.log(max_eig_B)))) if N == 0 else N\n",
    "\n",
    "    # iteration step\n",
    "    while N:\n",
    "        x = B@x_0+f\n",
    "        if q < 1 and q/(1-q)*F(x-x_0) < eps:\n",
    "            print('last iteration error: ',q/(1-q)*F(x-x_0))\n",
    "            break\n",
    "        x_0 = x\n",
    "        N -= 1\n",
    "\n",
    "    if test:\n",
    "        print(\"The speed of iteration is:{:.2f}\".format(-np.log(max_eig_B)))\n",
    "        \n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "id": "d519d2d9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "q = 0.41703337\n",
      "last iteration error:  4.200151608806934e-05\n",
      "The speed of iteration is:0.96\n",
      "[[-1.873]\n",
      " [ 0.694]\n",
      " [ 0.809]\n",
      " [-1.576]]\n"
     ]
    }
   ],
   "source": [
    "ans=Serial_iteration(A, b, N=0, eps=1e-4, iter_way='Gauss_Seidel',norm='infty', own_method=False, test=True)\n",
    "print(ans)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "id": "d93a18ce",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0009182348174605615"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Matrix_Solutions.m_infinite_norm(ans-x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "id": "954ce091",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.872],\n",
       "       [ 0.693],\n",
       "       [ 0.809],\n",
       "       [-1.577]])"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "id": "1da775df",
   "metadata": {},
   "outputs": [],
   "source": [
    "def SOR_iteration(A, b, eps=1e-4, x_0=None, omega=1, norm='infty', own_method=False, test=False):\n",
    "    \"\"\"Serial_iteration:Jacobi and Gauss_Seidel\n",
    "\n",
    "    Args:\n",
    "        A ([np.darray]): [A]\n",
    "        b ([np.darray]): [b]\n",
    "        eps ([float], optional): [error threshold]. Defaults to 1e-4.\n",
    "        x_0 ([np.darray], optional): [initial x for iteration]. Defaults to None.\n",
    "        omega (float, optional): [iteration coefficient]. Defaults to 1 means to Gauss_Seidel method.\n",
    "        norm (str, optional): [norm method to calculate fitting error,'infty','L1']. Defaults to 'infty'.\n",
    "        own_method (bool, optional): [use own method to calculate inv or norm, etc]. Defaults to False.\n",
    "        test (bool, optional): [show testing information]. Defaults to False.\n",
    "\n",
    "    Raises:\n",
    "        ValueError: [The spetrum of Iteration Matrix must be smaller than 1]\n",
    "\n",
    "    Returns:\n",
    "        [x]: [solution of Ax=b]\n",
    "    \n",
    "    Author: Junno\n",
    "        \n",
    "    Date: 2022-04-29\n",
    "    \"\"\" \n",
    "    assert A.shape[0] == b.shape[0]\n",
    "    assert 0<omega<2,'iteration coefficient must be during (0,2) to meet the convergence condition'\n",
    "    # use np method or your own\n",
    "    Inv = np.linalg.inv if not own_method else Matrix_Solutions.Primary_Row_Transformation_Inv\n",
    "    SVD = np.linalg.svd if not own_method else Matrix_Solutions.svd\n",
    "    # substract D,L,U from A\n",
    "    D, L, U = np.diag(np.diag(A)), np.tril(A, k=-1), np.triu(A, k=1)\n",
    "    # calculate B and f\n",
    "    D_=Inv((D+omega*L))\n",
    "    B =D_@((1-omega)*D-omega*U)\n",
    "    f=omega*D_@b\n",
    "    # initial x for iteration\n",
    "    if x_0 == None:\n",
    "        x_0 = np.zeros_like(b)\n",
    "    \n",
    "    # get eig of B\n",
    "    u, s, v = SVD(B)\n",
    "    # spectrum of B\n",
    "    eig_B = np.diag(s)\n",
    "    max_eig_B = np.max(eig_B)\n",
    "    # check convergence condition\n",
    "    if max_eig_B > 1:\n",
    "        print('The spetrum of Iteration Matrix must be smaller than 1')\n",
    "        raise ValueError\n",
    "        \n",
    "    # matrix norm\n",
    "    if norm == 'infty':\n",
    "        F = Matrix_Solutions.m_infinite_norm\n",
    "    elif norm == 'L1':\n",
    "        F = Matrix_Solutions.m_L1_norm\n",
    "\n",
    "    q = F(B)\n",
    "    if test:\n",
    "        print('q =', q)\n",
    "    \n",
    "    # iteration times\n",
    "    least_steps = int(np.ceil(np.log10(1/eps)*np.log(10) /\n",
    "            (-np.log(max_eig_B))))\n",
    "    N=0\n",
    "    # iteration step\n",
    "    while N<least_steps:\n",
    "        x = B@x_0+f\n",
    "        if q < 1 and q/(1-q)*F(x-x_0) < eps:\n",
    "            break\n",
    "        x_0 = x\n",
    "        N += 1\n",
    "\n",
    "    if test:\n",
    "        print(\"Iteration times:{}\".format(N))\n",
    "        print('last iteration error: ',q/(1-q)*F(x-x_0))\n",
    "        print(\"The speed of iteration is:{:.2f}\".format(-np.log(max_eig_B)))\n",
    "        \n",
    "    return x,N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "id": "e8ce08ec",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "q = 1.1293099524699077\n",
      "Iteration times:31\n",
      "last iteration error:  -0.0\n",
      "The speed of iteration is:0.31\n",
      "[[ 0.178]\n",
      " [-0.471]\n",
      " [-0.099]\n",
      " [-0.02 ]]\n"
     ]
    }
   ],
   "source": [
    "ans,N=SOR_iteration(A, b, eps=1e-4, omega=0.6,norm='infty', own_method=False, test=True)\n",
    "print(ans)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "id": "3dd97a30",
   "metadata": {},
   "outputs": [],
   "source": [
    "ws=np.arange(0.1,1.2,0.1)\n",
    "N_with_w=[]\n",
    "for w in ws:\n",
    "    ans,N=SOR_iteration(A, b, eps=1e-4, omega=w,norm='infty', own_method=False, test=False)\n",
    "    N_with_w.append(N)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "id": "9a4d991d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[182, 89, 58, 43, 35, 31, 28, 29, 32, 42, 68]"
      ]
     },
     "execution_count": 115,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N_with_w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "id": "40c09ea5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEXCAYAAABCjVgAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAeb0lEQVR4nO3de5wcZZ3v8c8XEkIgXM0wQu6EAJu4Am6HE3FXIoggi0Q5ejZkkYgoRw6rrLuEi+shBGQPctxzFDioEQPRFTAgShAEMSzga5fbhHsCSCRAhkAyyHW5hAR+54+qVJqhZ6Zmpqsr0/N9v17z6q6nquv5PT0z/et6nqp6FBGYmZkBbFF2AGZmtvlwUjAzs4yTgpmZZZwUzMws46RgZmYZJwUzM8s4KZiZWcZJwRpG0jJJ05utrv7oKU5JT0r6eD/2f5mkb9WqS9Jeku6T9Kqkr3Ve7mudNrA5KQwi1R8w/f2w6U1dG0XElIi4dSDXVW+d4yzy91LjPTkVuDUitouIC2osN0zRf4+Wn5OC9ZqkIWXHYHUxDljWzXIu/ntoLk4Kg5CknwJjgesk/aekU9Py3ST9QlKHpJXVXQjpN7nTJD0IvCZpiKTTJf0x7W5YLukzPey/+kjlzyTdKumltFvjyE4xPinpFEkPSnpZ0s8lbd2LtnQ+KpqT7us1ST+W1CrpN2nsv5O0U9U+u3sfTpP0TPq6xyQdXCOm4yRdV7W8QtKiquVVkvatEWfNtgD79vQ+VO17P0n3pvH9HNi6al11XbcAHwMuSuvqvLxnH/4euty+u99pN+3u3Lbt0vf+wE7lYySFpPd19b5YL0SEfwbJD/Ak8PHOz9PlLYClwJnAVsDuwBPAoVXb3w+MAYanZZ8Ddktf+zfAa8CutfZfXQYMBVYA30jrOgh4Fdir07Z3p/vfGXgE+EpP7eqmrXcCrcAoYC1wL7AfMAy4BZjb0/sA7AWsAnZLtx0PTKwRz+7AS+m+dgWeAp6pWvcisEUXv4day3nfh63Sur6evsefBdYD3+pi37cCX6q13Nu/h56276kttX6HNdp3FnBTF+v+E5he9v9YM/z4SME2mgq0RMTZEfFWRDwB/AiYWbXNBRGxKiLeAIiIqyJidUS8ExE/Bx4H9s9R1zRgBHBeWtctwK+Bozttd0G6/xeA64B9+9G+CyNiTUQ8A/weuCsi7ouIdcAvSRIEdP8+vE2SRCZLGhoRT0bEHztXlL7m1TTeA4GbgGck7Z0u/z4i3ulF7Hnfh2kkyeC7EbE+Iq4G7ulFPdV6+/eQZ/vetOVdJG0JnAhcki63SNq9apMNJMnJ+sl9gbbROGA3SS9VlW1J8gG60arqF0g6FvgHkm/MkHzQj8xR127Aqk4fjE+RfIuv9lzV89fT1/XVmqrnb9RYHpE+7/J9iIgVkv6e5BvrFEk3Af8QEatr1HcbMB3YI33+EklC+HC63Bt534fdSI5Iqm99/FQv69qot38PebaHvv9OPwDsQpJgIfm7E3C6pOHAdiRHgNZPTgqDV+d7pq8CVkbEpDyvkTSO5JvgwcAdEfG2pPtJ/lFr7b/aamCMpC2qEsNY4A+9iL9mXHXQ7fsQEZcDl0vaHvgh8G3g8zU2vQ34FDAB+GeSpPC3JEnhom7q709bngVGSVJVYhgLvOdoJode/T3k3L47PbV7FPBiRLySLh8G/DR9fiBJl9x9fazbqrj7aPBaQ9Lvu9HdwCvp4OFwSVtK+oCkqV28fluSf+QOSAZXSb7NdbX/aneRjD+cKmmoknPnPwVcWae29EeX74OS8/gPkjQMeJPkCOPtLvZzG8nA7fCIaCf5xnwY8D66//DqT1vuIOlG+Vo68HsU+brzaunt30Nvt++sp3a/AGwvaYKko0nGLSZL2pHkyO27veySsy44KQxe/wv4ppKzf06JiLdJPpj3BVYCz5P03+5Q68URsRz4F5IPojXAnwP/3tX+O732LeBI4JNpPRcDx0bEo/VoSx/3sTG27t6HYcB5adlzJN0Z3+hiP38gGfz8fbr8CsnA67+nddS9Len7ehTwBZJvzn8DXNObfVTtq7d/D73avoae2n0PyZeG+4HjSf5+DiAZx7qL5IjN6kDv7n40M7PBzEcKZmaWcVIwM7OMk4KZmWWcFMzMLDOgr1MYOXJkjB8/vuwwzMwGlKVLlz4fES211g3opDB+/Hja2trKDsPMbECR1OWV7u4+MjOzjJOCmZllnBTMzCzjpGBmZhknBTMzywzepPDyyzBlSvJoZmbAYE4K118Py5fDDTeUHYmZ2WZj8CWFWbNgxAiYPTtZPvbYZHnWrHLjMjPbDBSWFCQtkLRW0sNVZftKulPS/ZLaJO1fte4MSSskPSbp0KLi4uyzYexYGDo0WR46FMaNg3POKaxKM7OBosgjhctIZpqqdj4wLyL2Bc5Ml5E0mWSC7ynpay5OJ+quvz32SBLD+vWw7bbJ47x5MHFiIdWZmQ0khSWFiLidZAq9dxUD26fPdyCZqxdgBnBlRKyLiJXACvo+jWDPFi1KEsK8ecnjVVcVVpWZ2UDS6Hsf/T1wk6TvkCSkA9LyUcCdVdu1p2XvIekE4ASAsWPH9i2KOXPgwguhtRWOOQZWrerbfszMmkyjB5pPBL4eEWOArwM/TstVY9ua84RGxPyIqEREpaWl5k3+ejZ1apIQIHmsVPq2HzOzJtPopDCbTROJX8WmLqJ2YEzVdqPZ1LVkZmYN0uiksBo4MH1+EPB4+nwxMFPSMEkTgEnA3Q2Ozcxs0CtsTEHSFcB0YKSkdmAu8GXge5KGAG+Sjg1ExDJJi4DlwAbgpIh4u6jYzMystsKSQkQc3cWqv+hi+3OBc4uKx8zMejb4rmg2M7MuOSmYmVnGScHMzDJOCmZmlnFSMDOzjJOCmZllnBTMzCzjpGBmZhknBTMzyzgpmJlZxknBzMwyTgpmZpZxUjAzs4yTgpmZZZwUzMws46RgZmaZwpKCpAWS1kp6uFP5VyU9JmmZpPOrys+QtCJdd2hRcZmZWdcKm3kNuAy4CPjJxgJJHwNmAB+MiHWSdknLJwMzgSnAbsDvJO3pKTnNzBqrsCOFiLgdeKFT8YnAeRGxLt1mbVo+A7gyItZFxEpgBbB/UbGZmVltjR5T2BP4K0l3SbpN0tS0fBSwqmq79rTsPSSdIKlNUltHR0fB4ZqZDS6NTgpDgJ2AacAcYJEkAaqxbdTaQUTMj4hKRFRaWlqKi9TMbBBqdFJoB66JxN3AO8DItHxM1XajgdUNjs3MbNBrdFL4FXAQgKQ9ga2A54HFwExJwyRNACYBdzc4NjOzQa+ws48kXQFMB0ZKagfmAguABelpqm8BsyMigGWSFgHLgQ3AST7zyMys8ZR8Jg9MlUol2trayg7DzGxAkbQ0Iiq11vmKZjMzyzgpmJlZxknBzMwyTgpmZpZxUjAzs4yTgpmZZZwUzMws46RgZmYZJwUzM8s4KZiZWcZJwczMMk4KZmaWcVIwM7OMk4KZmWWcFMzMLOOkYGZmGScFMzPL9JgUJJ0vaXtJQyUtkfS8pGNyvG6BpLXp1Jud150iKSSNrCo7Q9IKSY9JOrT3TTEzs/7Kc6TwiYh4BTgCaAf2BObkeN1lwGGdCyWNAQ4Bnq4qmwzMBKakr7lY0pY56jAzszrKkxSGpo+HA1dExAt5dhwRtwO1tv2/wKlA9eTQM4ArI2JdRKwEVgD756nHzMzqJ09SuE7So0AFWCKpBXizL5VJOhJ4JiIe6LRqFLCqark9Lau1jxMktUlq6+jo6EsYZmbWhR6TQkScDnwYqETEeuB1km/2vSJpG+CfgDNrra5VdRfxzI+ISkRUWlpaehuGmZl1I89A8zbAScD306LdSI4aemsiMAF4QNKTwGjgXknvJzkyGFO17WhgdR/qMDOzfsjTfXQp8BZwQLrcDnyrtxVFxEMRsUtEjI+I8el+PhQRzwGLgZmShkmaAEwC7u5tHWZm1j95ksLEiDgfWA8QEW9Qu7vnXSRdAdwB7CWpXdLxXW0bEcuARcBy4EbgpIh4O0dsZmZWR0NybPOWpOGkffySJgLrenpRRBzdw/rxnZbPBc7NEY+ZmRUkT1KYS/LtfYyknwEfAb5QZFBmZlaOHpNCRNws6V5gGkm30ckR8XzhkZmZWcPlvffRKGBLYCvgo5KOKi4kMzMrS49HCpIWAB8ElgHvpMUBXFNgXGZmVoI8YwrTImJy4ZGYmVnp8nQf3ZHesM7MzJpcniOFhSSJ4TmSU1EFRER8sNDIzMys4fIkhQXA54GH2DSmYGZmTShPUng6IhYXHomZmZUuT1J4VNLlwHVUXckcET77yMysyeRJCsNJksEnqsp8SqqZWRPKc0XzcY0IxMzMytdlUpB0akScL+lCakx4ExFfKzQyMzNruO6OFB5JH9saEYiZmZWvy6QQEdelT1+PiKuq10n6XKFRmZlZKfJc0XxGzjIzMxvguhtT+CRwODBK0gVVq7YHNhQdmJmZNV53RwqrScYT3gSWVv0sBg7taceSFkhaK+nhqrL/LelRSQ9K+qWkHavWnSFphaTHJPW4fzMzq7/uxhQeAB6QdHlErO/Dvi8DLgJ+UlV2M3BGRGyQ9G2SbqjT0hvuzQSmALsBv5O0p+dpNjNrrB7HFPqYEIiI24EXOpX9NiI2dj3dCYxOn88AroyIdRGxElgB7N+Xes3MrO/yzrxWhC8Cv0mfjwJWVa1rT8veQ9IJktoktXV0dBQcopnZ4FJKUpD0TySD1T/bWFRjs/dcMAcQEfMjohIRlZaWlqJCNDMblPJMx7knMAcYV719RBzUlwolzQaOAA6OiI0f/O3AmKrNRpMMdJuZWQPluSHeVcAPgB8B/Rr4lXQYcBpwYES8XrVqMXC5pP9DMtA8Cbi7P3WZmVnv5UkKGyLi+73dsaQrgOnASEntwFySs42GATdLArgzIr4SEcskLQKWk3QrneQzj8zMGk+benC62EA6C1gL/JJ3z6fwQlevaZRKpRJtbb41k5lZb0haGhGVWuvyHCnMTh/nVJUFsHt/AzMzs81LnvkUJjQiEDMzK1+es4+GAicCH02LbgV+2NeL2szMbPOVp/vo+8BQ4OJ0+fNp2ZeKCsrMzMqRJylMjYh9qpZvkfRAUQGZmVl58lzR/LakiRsXJO1OP69XMDOzzVOeI4U5wL9JeoLkdhTjgOMKjcrMzEqR5+yjJZImAXuRJIVHI2JdDy8zM7MBqLuZ1w6KiFskHdVp1URJRMQ1BcdmZmYN1t2RwoHALcCnaqwLwEnBzKzJdDfz2tz06dnpxDcZSb6gzcysCeU5++gXNcqurncgZmZWvu7GFPYmmTN5h07jCtsDWxcdmJmZNV53Ywp7kUyGsyPvHld4FfhygTGZmVlJuhtTuBa4VtKHI+KOBsZkZmYlyXPx2n2STiLpSsq6jSLii4VFZWZmpcgz0PxT4P3AocBtJPMnv9rTiyQtkLRW0sNVZTtLulnS4+njTlXrzpC0QtJjkg7tfVPMzKy/8iSFPSLifwKvRcRC4K+BP8/xusuAwzqVnQ4siYhJwJJ0GUmTgZkkRyOHARdL2jJXC8zMrG7yJIWN8ya8JOkDwA7A+J5eFBG3A52n7JwBLEyfLwQ+XVV+ZUSsS6+JWAHsnyM2MzOrozxJYX7azfNNYDGwHPh2H+trjYhnAdLHXdLyUcCqqu3a07L3kHSCpDZJbR0dHX0Mw8zMaul2oFnSFsArEfEicDvFzcusGmVRa8OImA/MB6hUKjW3MTOzvun2SCEi3gH+ro71rZG0K0D6uDYtbwfGVG03Glhdx3rNzCyHPN1HN0s6RdKY9OyhnSXt3Mf6FgOz0+ezgWurymdKGpbeV2kScHcf6zAzsz7Kc53CxusRTqoqC3roSpJ0BTAdGCmpHZgLnAcsknQ88DTwOYCIWCZpEcl4xQbgpIjw7G5mZg2miIHbLV+pVKKtra3sMMzMBhRJSyOiUmtdj91HkraR9E1J89PlSZKOqHeQZmZWvjxjCpcCbwEHpMvtwLcKi8jMzEqTJylMjIjzSS9ii4g3qH0KqZmZDXB5ksJbkoaTXjcgaSKwrtCozMysFHnOPjoLuBEYI+lnwEeA44oMyszMytFjUoiI30paCkwj6TY6OSKeLzwyMzNruDxnHy2JiD9FxPUR8euIeF7SkkYE15RefhmmTEkezcw2M10mBUlbp1cuj5S0U9XVzOOB3RoWYbO5/npYvhxuuKHsSMzM3qO7I4X/DiwF9gbuTZ8vJbk1xf8rPrQmM2sWjBgBs9O7fBx7bLI8a1a5cZmZVekyKUTE9yJiAnBKREyo+tknIi5qYIzN4eyzYexYGDo0WR46FMaNg3POKTcuM7MqXQ40SzooIm4BnpF0VOf1EXFNoZE1mz32SBLD0UfDttvCunUwbx5MnFh2ZGZmme66jw5MHz9V48e3ueiLRYuShDBvXvJ41VVlR2Rm9i6+IV4j3XNP0oXU2gpr1sCqVVCpeU8qM7PCdHdDvDwXr1m9TJ266Xlra/JjZrYZyXObCzMzGyTyXLw2LE+ZmZkNfHmOFO7IWWZmZgNcd6ekvh8YBQyXtB+bbpe9PbBNfyqV9HXgSyR3Xn2I5AZ72wA/B8YDTwL/LSJe7E89ZmbWO90NNB8KfAEYDfwLm5LCq8A3+lqhpFHA14DJEfFGOjfzTGAysCQizpN0OnA6cFpf6zEzs97rMilExEJgoaT/GhG/KKDe4ZLWkxwhrAbOAKan6xcCt+KkYGbWUHnGFEZL2l6JSyTdK+kTfa0wIp4BvgM8DTwLvBwRvwVaI+LZdJtngV1qvV7SCZLaJLV1dHT0NQwzM6shT1L4YkS8AnyC5IP6OOC8vlYoaSdgBjCB5G6r20o6Ju/rI2J+RFQiotLS0tLXMMzMrIY8SWHjWMLhwKUR8QD9m6P548DKiOiIiPXANcABwBpJuwKkj2v7UYeZmfVBnqSwVNJvSZLCTZK2A97pR51PA9MkbSNJwMHAI8BiIL2vNLNJbtFtZmYNlOc2F8cD+wJPRMTrkt5HP+Zojoi7JF1NMkfDBuA+YD4wAlgk6XiSxPG5vtZhZmZ9k2eO5nckrQT2lLR1PSqNiLnA3E7F60iOGszMrCQ9JgVJXwJOJrle4X5gGskVzQcVGpmZmTVcnjGFk4GpwFMR8TFgP8DngpqZNaE8SeHNiHgTkhvhRcSjwF7FhmVmZmXIM9DcLmlH4FfAzZJeJLkC2czMmkyegebPpE/PkvRvwA7AjYVGZWZmpcg185qkvwQmRcSlklpI7p66stDIzMys4fJMsjOX5MZ0Z6RFQ4F/LTIoMzMrR56B5s8ARwKvAUTEamC7IoMyM7NuvPwyTJmSPNZZnqTwVkQEyYQ4SNq27lGYmVl+118Py5fDDTfUfdd5ksIiST8EdpT0ZeB3wI/qHomZmXVv1iwYMQJmp7eJO/bYZHnWrLpVkefso+9IOgR4heT6hDMj4ua6RWBmZvmcfTbcfz88+SRs2ABDh8K4cXDOOXWrItfZR2kScCIwMyvTHnskieHoo2HbbWHdOpg3DyZOrFsVXXYfSXpV0is1fl6V9ErdIjAzs/wWLUoSwrx5yeNVV9V1993N0ewzjMzMNjdz5sCFF0JrKxxzDKxaVdfd5+o+MjOzzcTUqZuet7YmP3WU5+wjMzMbJJwUzMwsU0pSkLSjpKslPSrpEUkflrSzpJslPZ4+7lRGbGZmg1lZRwrfA26MiL2BfYBHgNOBJRExCViSLpuZWQM1PClI2h74KPBjgIh4KyJeAmYAC9PNFgKfbnRsTa3Ae6WYWfMo40hhd5LpPC+VdJ+kS9L7KbVGxLMA6eMutV4s6QRJbZLaOjo8K2huBd4rxcyaRxlJYQjwIeD7EbEfyd1Xc3cVRcT8iKhERKWlpaWoGJtHA+6VYmbNo4yk0A60R8Rd6fLVJElijaRdAdLHtSXE1nzOPhvGjk3ukQKF3CvFzJpHw5NCRDwHrJK0V1p0MLAcWAykX2eZDVzb6Nia0sZ7paxfn1wSv3593e+VYmbNo6yzj74K/EzSg8C+wD8D5wGHSHocOCRdtnoo+F4pZtY8lMyfMzBVKpVoa2srO4zN3z33JF1Ira2wZk1yr5RKpeyozKwkkpZGRM0PAd/7aDAo+F4pZtY8fJsLMzPLOCmYmVnGScHMzDJOCmZmlnFSMDOzjJOCmZllnBTMzCzjpGBmZhknBTMzyzgpmJlZxknBzMwyTgpmZpZxUrDieX5oswHDScGK5/mhrRk16ZcdJwUrjueHtmbWpF92SksKkraUdJ+kX6fLO0u6WdLj6eNOZcVmdeL5oa0ZNfmXnTKPFE4GHqlaPh1YEhGTgCXpsg1knh/amlGTf9kpJSlIGg38NXBJVfEMYGH6fCHw6QaHZUXw/NDWbJr8y05ZRwrfBU4F3qkqa42IZwHSx11KiMvqbc4ceOwx+Md/TB7nzCk7IrP+a+IvOw2fo1nSEcDaiFgqaXofXn8CcALA2LFj6xuc1Z/nh7ZmNGcOXHhh8vd8zDGwalXZEdVNw5MC8BHgSEmHA1sD20v6V2CNpF0j4llJuwJra704IuYD8wEqlUo0Kmgzs0wTf9lpePdRRJwREaMjYjwwE7glIo4BFgPpcD6zgWsbHZuZ2WC3OV2ncB5wiKTHgUPSZbO+a9KLi8yKVGpSiIhbI+KI9PmfIuLgiJiUPr5QZmzWBJr04iKzIm1ORwpm9dHkFxeZFclJwZpPk19cZFYkJwVrPk1+cZF14rGjunJSsOZU9sVF/qBqHI8d1ZWTgjWnsq+k9gdV8Tx2VAgnBWtOU6duuqCotRUqlcbU6w+qxvHYUSGcFMzqaTB/UDW6y8xjR4VwUjCrp7I/qMocyyijy6zssaMm5KRgVm9lflCV8cFcZpdZ2WNHTUgRA/eecpVKJdra2soOw+zd7rkn6UJqbYU1a5I7aBY9pjFrFixeDOvWwYYNMGQIDBsGRx4Jl19ebN0rViT1PPkkvPEGDB8OEyYk8bgrZ7MkaWlE1Pyj9JGCWb2VMchd5lhG2V1mVldOCmbNoOwPZvftNw0nBbNmUeYHs/v2m4bHFMyaRRljGTYgdTemUMbMa2ZWhCaeDcwax91HZmaWcVIwM7OMk4KZmWWcFMzMLOOkYGZmmQF9SqqkDuCpsuPog5HA82UH0WBu8+Aw2No8UNs7LiJaaq0Y0ElhoJLU1tU5ws3KbR4cBlubm7G97j4yM7OMk4KZmWWcFMoxv+wASuA2Dw6Drc1N116PKZiZWcZHCmZmlnFSMDOzjJNCgSQdJukxSSsknV5j/d9KejD9+Q9J+5QRZ7301N6q7aZKelvSZxsZXxHytFnSdEn3S1om6bZGx1hvOf6ud5B0naQH0jYfV0ac9SRpgaS1kh7uYr0kXZC+Jw9K+lCjY6ybiPBPAT/AlsAfgd2BrYAHgMmdtjkA2Cl9/kngrrLjLrK9VdvdAtwAfLbsuBvwO94RWA6MTZd3KTvuBrT5G8C30+ctwAvAVmXH3s92fxT4EPBwF+sPB34DCJg2kP+XfaRQnP2BFRHxRES8BVwJzKjeICL+IyJeTBfvBEY3OMZ66rG9qa8CvwDWNjK4guRp8yzgmoh4GiAiBnq787Q5gO0kCRhBkhQ2NDbM+oqI20na0ZUZwE8icSewo6RdGxNdfTkpFGcUsKpquT0t68rxJN80Bqoe2ytpFPAZ4AcNjKtIeX7HewI7SbpV0lJJxzYsumLkafNFwJ8Bq4GHgJMj4p3GhFea3v6/b7Y881pxVKOs5vm/kj5GkhT+stCIipWnvd8FTouIt5MvkQNenjYPAf4COBgYDtwh6c6I+EPRwRUkT5sPBe4HDgImAjdL+n1EvFJwbGXK/f++uXNSKE47MKZqeTTJN6d3kfRB4BLgkxHxpwbFVoQ87a0AV6YJYSRwuKQNEfGrhkRYf3na3A48HxGvAa9Juh3YBxioSSFPm48Dzouks32FpJXA3sDdjQmxFLn+3wcCdx8V5x5gkqQJkrYCZgKLqzeQNBa4Bvj8AP7muFGP7Y2ICRExPiLGA1cD/2MAJwTI0WbgWuCvJA2RtA3wX4BHGhxnPeVp89MkR0ZIagX2Ap5oaJSNtxg4Nj0LaRrwckQ8W3ZQfeEjhYJExAZJfwfcRHLGxoKIWCbpK+n6HwBnAu8DLk6/PW+IAXrHxZztbSp52hwRj0i6EXgQeAe4JCJqntY4EOT8PZ8DXCbpIZJuldMiYiDeXjoj6QpgOjBSUjswFxgKWZtvIDkDaQXwOsnR0oDk21yYmVnG3UdmZpZxUjAzs4yTgpmZZZwUzMws46RgZmYZJwUzM8s4KZiZWcZJwayOJB0p6epOZSdKuqCsmMx6w0nBrL7OBc7qVPZHYHLjQzHrPScFszpJZ87bIiIeljRO0onpqqEM0Dtm2uDjpGBWP/sCS9PnhwCT0ueTSWYoM9vsOSmY1c8WwAhJWwJHkcw+Nhz4AnB5mYGZ5eWkYFY/N5DMXXw/yexyU4A2YH5E3FtiXGa5+S6pZmaW8ZGCmZllnBTMzCzjpGBmZhknBTMzyzgpmJlZxknBzMwyTgpmZpb5/ye6TboyaykYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplot(1,1,1)\n",
    "plt.scatter(ws,N_with_w,c='r',marker='*')\n",
    "plt.xlabel('$\\omega$')\n",
    "plt.ylabel('least iteration times')\n",
    "plt.title('Iteration times with different $\\omega$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "id": "48d879b6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.6472612633470305\n"
     ]
    }
   ],
   "source": [
    "D=np.diag(np.diag(A))\n",
    "U=np.triu(A,k=1)\n",
    "L=np.tril(A,k=-1)\n",
    "B = -np.linalg.inv(D)@(L+U)\n",
    "u, s, v = np.linalg.svd(B)\n",
    "wopt=2/(1+np.sqrt(1-np.max(np.diag(s))**2))\n",
    "print(wopt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82d6b36d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
